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ABSTRACT 

We describe new implementations of leptonic and hadronic models for 
the broadband emission from relativistic jets in AGN in a temporary steady 
state. For the leptonic model, a temporary equilibrium between particle injec- 



i-H I tion/acceleration, radiative cooling, and escape from a spherical emission region is 

T I evaluated, and the self-consistent radiative output is calculated. For the hadronic 

Q^l model, a temporary equilibrium between particle injection/acceleration, radiative 

Q I and adiabatic cooling, and escape is evaluated for both primary electrons and 

protons. A new, semi-analytical method to evaluate the radiative output from 

c^ I cascades initiated by internal 77 pair production is presented. We use our codes 

to fit snap-shot spectral energy distributions of a representative set of Fermi- 

^ I LAT detected blazars. We find that the leptonic model provides acceptable fits 

Y^ [ to the SEDs of almost all blazars with parameters close to equipartition between 

\^ I the magnetic field and the relativistic electron population. However, the hard 

^ I 7-ray spectrum of AO 0235-1-164, in contrast to the very steep IR-optical-UV 

^ I continuum, poses a severe problem for the leptonic model. If charge neutrality 

rn ' in leptonic models is provided by cold protons, the kinetic energy carried by the 

jet should be dominated by protons. We find satisfactory representations of the 

snapshot SEDs of most blazars in our sample with the hadronic model presented 

here. However, in the case of two quasars the characteristic break at a few GeV 

energies can not be well modelled. All of our hadronic model fits require powers 

in relativistic protons in the range Lp ~ 10^^ - 10^^ erg s~^. 
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Introduction 



Blazars are a class of radio-loud active galactic nuclei (AGNs) comprised of Flat- 
Spectrum Radio Quasars (FSRQs) and BL Lac objects. Their spectral energy distributions 
(SEDs) are characterized by non-thermal continuum spectra with a broad low-frequency com- 
ponent in the radio - UV or X-ray frequency range and a high-frequency component from 
X-rays to 7-rays, and they often exhibit substantial variability across the electromagnetic 
spectrum, in extreme cases on time scales down to just a few minutes. Blazars are sub-divided 
based on the location of their synchrotron peak: Low-Synchrotron-Peaked (LSP) blazars, 
consisting of FSRQs and Low-frequency-peaked BL Lac objects (LBLs) have z/gy < 10^^ Hz; 
Intermediate-Synchrotron-Peaked (ISP) blazars (which are exclusively BL Lac objects and 
hence also termed IBLs for Intermediate BL Lac Objects) have 10^^ Hz < i/gy < 10^^ Hz; 
High-Synchrotron-Peaked (HSP) blazars (also exclusively BL Lac objects, termed HBLs for 
High-peaked BL Lac objects) have Us > 10^^ Hz. 

The extreme inferred isotropic luminosities, combined with the rapid high-energy vari- 
ability, provide convincing evidence that the high-energy emission from blazars originates in 
relat ivistic jets closely aligned with our line of sight (for a review of those arguments, see, 
e.g.. ISchlickeiserlll996l ). The simplest (and often sufficient) assumption is that the emission 
is produced in an approximately spherical region, propagating along the jet with a speed 
(3rc, corresponding to a bulk Lorentz factor F. If the jet forms an angle 6'obs with respect 
to our line of sight, this results in Doppler boosting characterized by the Doppler factor 
D = (T [1 — (3t cos 6'obs])" • The observed bolometric flux will be enhanced by a factor D"^, 
while photon energies are blue-shifted by a factor D, and variability time scales will be short- 
ened by a factor D~^ (for a review of relativistic effects in the jets of AGN, see iBottcher 
2012h . 



It is generally accepted that the low-frequency (radio through UV or X-ray) emission 
from blazars is synchrotron emission from relativistic electrons in the jet. For the origin 
of the high-energy (X-ray through 7-ray) emission, two fundamentally different approaches 
have been proposed, generally referred to as leptonic and hadronic models. 

In leptonic models, the radiative output throughout the electromagnetic spectrum is as- 
sumed to be dominated by leptons (electrons and possibly positrons), while any protons that 
are likely present in the outflow, are not accelerated to sufficiently high energies to contribute 
significantly to the radiative output. The high-energy emission is then most plausibly ex- 
plained by Compton scattering of low-ener gy photons by the same electrons producing the 



synchrotron emission at lower frequencies ( Maraschi et al.l 1992 ■ Bloom fc Marscherl Il996 



Dermer fc Schlickeiserl Il993l : ISikora et al.lll994j : iBlazejowski et al.l 120001 ). In hadronic mod- 
els, both primary electrons and protons are accelerated to ultrarelativistic energies, with 
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protons exceeding the threshold for p7 photo-pion production on the soft photon field in the 
emission region. While the low-frequency emission is still dominated by synchrotron emission 
from primary electrons, the high-energy emission is dominated by proton synchrotron emis- 
sion, tt" decay photons, synchrotron and Compton emission from secondary decay products 
of charged pions, and the output from pair c ascades initiated by these high - energy emissions 



intrinsically absorbed by 77 pair production (JMannheim fc Biermanrull992l : lAharoniarul2000 



Miicke fc Protheroell200ll : iMiicke et al.ll2003r). For a gene ral overview of the features of both 
leptonic and hadronic blazar models, see iBottcherl ( 2010l ). 



Steady-state leptonic models have met with gr eat success in modeling the spectral energy 



Ghiselhni et al.lll998l: ICelotti fc Ghisellini 



Ghisellini et al.ll2010f ). 



distributions (SEDs) of all classes of blazars (e.g. 
20081 ). including substantial samples of Ferm2-detected blazars (e.g. 
Much progress has also been made to explore th e variability features predicted in time- 
dependent implementations of leptonic models (e.g..lMastichiadis fc Kirklll997l:lKusunose et al 



2OOOI: iLi fc KusunosdbOOd: iBottcher fc Chiang 



2002 



Sokolov et al.l 



2004; 



Graff et al. 



2008 



Bottcher fc Dermeii I2OIOI : IJoshi fc Bottcherl I2OII ). However, the very fast (time scal es of a 
few minutes) variability of some TeV blazars ( Albert et al.ll2007l : lAharonian et al.ll2007l ) poses 
severe problems for single- zone mode ls of blazars due to the required extremely high bulk 
Lorentz factors (JBegelman et al.ll2008[ ). Even withou t considering variab ility, the SED of the 
quasar 3C 279, detected in VHE 7-rays by MAGIC ( lAlbert et al.ll2008l ). pos es problems for 



one- z one leptonic models, and is more easily represented by a hadronic model (IBottcher et al. 
20091 ). In order to remedy some of the problems, several variations of multi-zo n e mod els 



have been proposed, includin g the spine-sheath model of 



the decelera ting-iet model of iGeorganopoulos fc KazanasI (20031). as well as internal-shock 



models (e.g..lMarscher fc Geaiill985l:ISpada et al. 



Joshi &: BottcheJbOllI : IBottcher fc Dermei]l2010[ ). 



2001 



^avecchio fc Ghiselhnil (|2008[ ) or 



Sokolov et al.ll2004J : lGraff et al.ll2008 



The goal of this work is to provide tools for the modeling of blazar spectra including 
data from Fermi-LAT , which have been accumulated over one year of Fermi operations. 
Any short-term variability in the SEDs of the blazars under consideration, has therefore 
been averaged out over this integration time. Therefore, for the purpose of this study, we 
will use time-independent models, which represent an appropriate time average of the rapidly 
variable broadband emission. 

The accurate evaluation of the radiative outpu t of hadronic models is best achieved by 



Monte-Carlo simulations (e.g., iMiicke et al.ll2000al ) due to the comphcated energy depen- 



dence of the cross sections (in particular, for p7 interactions) involved. As those are quite 
time-consuming, hadronic models have so far received much less attention in the literature 
than leptonic ones, and in particular, time dependent implementations of hadronic models 
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are still in their early development stages. A comparison between the results of fitting a 
large, statistically representative sample of blazars with both leptonic and hadronic models, 
considered to comparable degree of detail, has therefore never been performed. This is the 
main purpose of the code developments and modeling efforts presented in this paper. 

In ^we describe recent developments to our existing leptonic radiation transfer code, in 
particular its application to quasi-stationary situations, and the implementation of arbitrary 
external radiation fields as sources for Compton scattering. In §21 we will describe a new, 
semi-analytical implementation of the stationary hadronic model and test it against r esults 
of detailed Monte-Carlo simulations based on the SOPHIA code (JMiicke et al.ll2000al ). We 
present the results of our comparative leptonic and hadronic modeling of the SEDs of a 
representative sample of Fermi-LAT detected blazars of different sub-classes in ^ We 
summarize and discuss our results in §3 Throughout our paper, we convert redshifts to 
luminosity distances using a ACDM cosmology with Hq = 70 km s^^ Mpc"^, Qm = 0.3 and 
^A = 0.7. 



2. Stationary leptonic blazar model 

The leptonic jet radi ation tra nsfer model used for this work is based on the work of 
Bottcher fc Chiang! (120021 ) (see also lBottcher. Mause &: Schlickeiserlll997l : iBottcher &: Bloom 
2000). It is a homogeneous one- zone model in which a population of ultrarelativistic electrons 
(or positrons) is injected with a power-law distribution Qdle) = Qo le'^" H^^je', 7e,i, 7e,2) into 
a spherical emission region of co-moving radius R. Here, H{x; a, b) is the Heaviside function 
defined asH = liia<x<b and H = otherwise. The emission region moves with 
constant relativistic speed (S^c (corresponding to bulk Lorentz factor F) along the jet. 

The electron distribution cools due to synchrotron and Compton emission. Synchrotron 
emission is determined through a tangled magnetic field of co- moving strength B. For 
Compton scattering, the synchrotron radiation field (SSC = synchrotron self-Compton) and 
various external radiation fields are taken into account: (a) Direct accretion disk emission (see 



Bottcher. Mause fc Schlickeiseii Il997l . for details), (b) accretio n disk emission re-processed 
by the Broad Line Region (BLR; see iBottcher fc Bloomll2000l . for details), (c) an isotropic 
(in the rest frame of the AGN) external radiation field of arbitrary spectral shape. For 
the latter, we have developed a pre-processing routine which produces a spectrum file of 
the external radiation field in the AGN rest frame. The blazar code applies the proper 
relativistic Lorentz transformations into the blob rest frame to use this radiation field as 
sources for external Compton scattering. Such a description is appropriate for radiation 
fields for which the angular characteristics of the radiation field in the co-moving frame are 
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dominated by relativistic aberration rather than any intrinsic anisotropy. Quantitatively, 
the characteristic scale of angular variations A6'ext of the external radiation field should be 
A^ext ^ l/f- This is typically the case for (1) (line-dominated) emission from the BLR as 
long as the emission region is located within the BLR; (2) infrared emission from a large- 
scale dusty torus around the central engine; (3) the Cosmic Microwave Background. The 
external Compton emissivity is evaluated using the head-on approximation for the Compton 
cross section (based on an angle integration of the full Klein-Nishina cross section; see, 



e.g., iDermer fc Menoru |2009| ). The electron cooling rates due to Compton scattering are 
calculated using the full Klein - Nishin a cross section, adopting the analytical solution of 
Bottcher. Mause fc Schlickeiserl ( 119971 ). Particles may escape the emission region on a time 
scale tesc = Vesc R/c, which we parameterize with an escape time scale parameter ?7esc > 1- 

In order to fit SEDs of blazars in the absence of detailed spectral variability information, 
it is appropriate (and computationally more efficient than detailed time-dependent radiation 
transfer simulations) to apply a stationary radiation transfer model. For the stationary model 
used in this work, our code finds an equilibrium between the relativistic particle injection 
mentioned above, radiative cooling, and escape. For this purpose, a critical Lorentz factor 
7c is determined, for which the radiative cooling time scale equals the escape time scale, 
7c/|7(7c)| = tesc- The shape of the resulting equilibrium particle distribution will then 
depend on whether •jc < 7e,i (the fast cooling regime) or 7c > 7e,i (the slow cooling regime). 
In the fast cooling regime. 



7e 
7<:,1 



^e^'ile) = no < 







for 7c < 7e < 7e,l 



^ ) for 7c,l < 7e < 7e,2 



else 



while in the slow cooling regime. 



nTi7e) = no< 







-fe+i) 



for 7e,l < 7e < 7c 



for 7c < 7e < 7e,2 



else 



(2) 



Since the radiative cooling depends on the self-consistent radiation field, including syn- 
chrotron emission from the present relativistic electron population, an iterative scheme is 
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applied: The code starts out with an equihbrium based on only synchrotron and external- 
Compton coohng. This distribution is then used to calculate the synchrotron radiation 
energy density, which is added to the total radiation energy density to re-calculate 7c and 
the resulting equilibrium particle distribution. The process is iterated until convergence is 
achieved. 

The code then evaluates the resulting kinetic power in relativistic electrons in the AGN 
frame, 

oo 
Le = IT R^ r^ /3r C TTleC^ / d'ye He {'Je) 7e (3) 

1 

and compares it to the power carried in the magnetic field (Poynting flux), 

LB = nR'T'(3rc^ (4) 

We deflne the equipartition parameter as the ratio of the two, i.e., ese = Ls/Le- We also 
evaluate the kinetic luminosity in cold protons that are expected to be present if charge 
neutrality is provided by one cold proton per electron (i.e., no pairs): 

Lp = 71 R'^T'^ (3rcnempC^ (5) 

The presence of pairs in addition to an electron-proton plasma in the jet will lower the total 
kinetic luminosity of the jet, Lkin = L^ + Lp. 

The solid line in Figure [1] shows a generic model blazar SE P with an EC component 



calculated with the full BLR re-processing calculation described in lBottcher fc Bloom! (J200d ) 
(dot-dashed line) for a blazar located at a redshift of 2; = 0.3. For this calculation, we assumed 
a Shakura-Sunyaev accretion disk with a luminosity of L^ = 10'^^ erg s~^ with a multi-color 
blackbody spectrum (the dotted line in Figured]) peaking at u^ ~ 3.9 x 10^^ Hz in the AGN 
rest frame. A spherical BLR with re-processing fraction t^blr = 0.04 is placed at a distance 
-Rblr = 0.2 pc from the accretion disk. The average radiation energy density inside the BLR 
is commonly estimated as 

ublr = -. — ^2 ~ 2.9 X 10 erg cm (6) 

In an emission region located close to the inner edge of the BLR, an electron distribution 
is injected with a power-law with index q = 2.5 between 71 = 10^ and 72 = 5 x 10^. The 
equilibrium electron population is in the fast cooling regime. 
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Fig. 1. — Generic blazar SEDs with EC (BLR) calculated with the full BLR geometry (solid 
SED; EC [BLR] component: dot-dashed line). The dotted line shows the disk spectrum. 
The dashed line shows the EC component calculated using an isotropic, thermal external 
radiation field with a temperature producing the same peak frequency as the disk, and the 
same radiation energy density as resulting from BLR reprocessing of the disk emission. 



For comparison, the dashed line in Figure [T] shows the (much faster) calculation of 
the EC component using an isotropic, thermal radiation field with a blackbody temperature 
T = 5000 K, reproducing a peak energy close to the one of accretion disk, and with an energy 
density of itext = 2.5 x 10~^ erg cm~^. The figure shows that the isotropic radiation field 
approximation with these parameters provides an excellent description of the EC component. 
This indicates that the approximation of Eq. ([H]) slightly over-estimates the energy density 
of the radiation field. In the comprehensive modeling effort described in §U we represent all 
BLR components as isotropic radiation fields, keeping in mind that a simple translation of 
disk and BLR parameters to the external radiation field energy density of Eq. ([6]) needs to 
be taken with caution. 

The model described above has already been used successfully to model a number of 
blazar SEDs. In particular, it has been used to interpret the SEDs of blazars detected 
in very high energies (VHE, E > 100 GeV) by t he Very Energetic Radiatio n Imaging Tele- 
scope Array System (VERITAS), e.g., W Comae jAcciari et al.ll2008l . l2009bl ). lES 0806+524 



flAcciari et alJl2009al) . PKS 1424+240 (lAcciari et alJl2010ah . RGB J0710+591 flAcciari et al. 
2010^, and 3C 66A jAbdo et al.l |201l|). It has also been used to model the SED of the 



high-redshift FSRQ PKS 0528+134 in quiescence flPalma et al.ll201l[ ). 



3. Stationary hadronic blazar model 



In order for protons to contribute significantly to t 



le radiative output in relativistic jets 



of bla zars either through proton synchrotron emission (|A 



laroman 



2000 



Miicke &: Protheroe 



20001 ) or photo-pion production (JMannheim &: Biermannlll992l ). they need to be accelerated 
to proton energies of typically Ep = Eiq 10^^ eV with Eiq > 1. Requiring that the Larmor 
radius of these protons be smaller than the size of the emission region, R = 10^^ R15 cm, 
one needs magnetic fields of -B > 30 Eiq R^^ G. The energy density of such large magnetic 
fields is likely to dominate over the energy density of external radiation fields in typical 
AGN environments, so that the synchrotron radiation field is likely to be the dominant 
target photon field for photo-pion production by ultrarelativistic protons. We will therefore 
restrict our considerations in thi s pape r to t he synchrotron-prot on blazar model as discussed 
in detail in iMiicke fc Protheroel feoOlJ ) and iMiicke et all J2003I ). 



For a comparison between leptonic and hadronic blazar models on the basis of a suffi- 
ciently large sample, we need a hadronic model implementation that allows for the efficient 
scanning through parameter space within a reasonable time frame. At the same time, the 
model needs to treat hadronic processes to a comparable degree of detail as the leptonic pro- 
cesses are being considered in the leptonic model described in the previous section. The most 
accurate and reliable method to evaluate the radiative output from photo-pair and photo- 
pion production and the cascades initiated by 77 absorption of ultra-high-energy (UHE; 
E ^ 1 TeV) photons produced through 7r° decay and synchro tron emission of pio ns and 



their decay products, is through Monte-Carlo simulations (e.g.. IMiicke et al.ll2000al ). How- 
ever, such simulations are quite time consuming, and a comprehensive modeling effort of a 
large sample of SEDs is currently infeasible with this method. 



In order to avoid time-consuming Monte-Carlo simulations, iKelner fc AharonianI (J2008[ ) 
have developed analytical expressions for the final decay products (electrons, positrons, pho- 
tons, neutrinos) from photo-pion production for isotropically distributed, mono-energetic 
photons and protons. These can be integrated over any proton and photon distribution to 
obtain the final production spectra of electrons, positrons, photons, and neutrinos. However, 
the 7r° decay photons as well as synchrotron photons from these first-generation electrons 
and positrons emerge in the UHE regime, at which the dense radiation fields in the emission 
regions of blazars are highly opaque to 77 pair production. It is therefore essential to include 
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the effects of UHE 7-ray induced pair cascades in hadronic models. Again, time-consuming 
Monte- Carlo/numerical simulations constitute the standard method for the evaluation of 
such cascades. 



3.1. Synchrotron-Supported Pair Cascades 

For the purpose of an efficient calculation of the cascades, we have developed a semi- 
analytical method that does not involve Monte-Carlo simulations. Let us assume that the 
injection rates of first-generation high- energy 7-rays, N^, and pair s, Qdl), are known, e.g.. 



from the analytical approximations of iKelner fc Aharoniaru (J2008[ ). Throughout this expo- 



sition, we use the dimensionless photon energy e = hv/rUeC . In the case of linear cascades 
the optical depth for 77 absorption, T^7(e) can be pre-calculated from the low-energy radi- 
ation field. Under these conditions, the spectrum of escaping (observable) photons can be 
calculated as 

where N^^ has contributions from the first-generation high-energy photon spectrum and 
synchrotron emission from secondaries, iV^^™ = A^° + A*"^^^. We use a synchrotron emissivity 
function for a single electron of the form jy oc z/^/^ g-^/^o ^jth ^q = 67^, where b = B/Bcnt 
and Ecrit = 4.4 x lO^^ G. This yields 



1 
with the normalization 

carB"^ 
° " evrmec^r (4/3) 64/3 ^^' 

The electron distribution, Ne{'~f) will be the solution to the isotropic Fokker-Planck 
equation in equilibrium {d{.)/dt = 0): 

^ (7 Ar,[7]) = Q,(7) + Nril) + NeilT''- (10) 

In the case under consideration here, the electron energy losses will be dominated by syn- 
chrotron losses, i.e., 
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carB'^ 2— 2 n^\ 

7 = -7 ^7 =-1^01- 11 

ovr rrieC^ 

We represent the escape term in Eq. (ITOl) through an energy-independent escape time scale 
^esc = VescR/c, SO that Nei^^Y^^ = " ^e(7)Aesc- ^e^ il) i^ Eq. [TU] is the rate of particle 
injection due to 77 absorption, to be evaluated self-consistently with the radiation field. In 
the 77 absorption of a high-energy photon of energy e, one of the produced particles will 
assume the major fraction, /^ of the photon energy. From comparison with Monte-Carlo 
simulations, we find that /^ = 0.9 yields good agreement with numerical solutions to the 
cascade problem. Hence, an electron/positron pair with energies 71 = /^ e and 72 = (1 — /-y) e 
is produced. Furthermore realizing that every photon not escaping (according to Eq. [7]) will 
produce an electron/positron pair, we can write the pair production rate as 

Nrii) = /ab.(ei) (iv° + n:i) + Us{^,) (iv° + n:i) (12) 

where ei = 7//^, 62 = 7/(1 - A) and 

/abs(e) = 1 '—-- (13) 

With this approximation, we find an implicit solution to Equation [TOl 

00 
Neil) = A /^7 (ge(7) + Nril) - ^1 (14) 

7 

The solution JHM is implicit in the sense that the particle spectrum A^e(7) occurs on both 
sides of the equation as N^'^ depends on the synchrotron emissivity calculated through 
Eq. [HI which requires knowledge of N^i'^), where pairs at energies of 71 = ^/'y/if-yb) and 
72 = a/7/([1 — f-y]b) provide the majority of the radiative output relevant for pair production 
at energy 7. However, for practical applications, we may use the fact that generally, 7, the 
argument on the l.h.s., is much smaller than 71^2- Therefore, Eq. [14] may be evaluated 
progressively, starting at the highest pair energies for which Qo(7) 7^ or A^^,^^ ¥" O5 and 
then using the solution for Ne{'y) for large 7 as one progresses towards lower values of 7. 

Once the equilibrium pair distribution Ne{'~f) is known, it can be used in Eq. [8] to 
evaluate the synchrotron emissivity and hence, using Eq. [7] the observable photon spectrum. 

Figure [2] illustrates the Monte-Carlo generated synchrotron + cascade spectra (solid 
lines) for the injection of monoenergetic electrons with energies 70 = 10^°, 10^^ and 10^^, 
respectively. The magnetic field is i? = 10 G, and the 77 optical depth as a function of 
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Fig. 2. — Comparison of the cascade emission from Monte-Carlo simulations (solid curves) 
and our semi-analytical description (dashed curves) in the case of the injection of mono- 
energetic electrons with energies given by the labels, in a -B = 10 G magnetic field. The blue 
curve shows the 77 opacity as a function of 7-ray photon frequency. 

photon energy is shown by the blue solid curve. The dashed lines illustrate the results of the 
analytic approximations developed here. The agreement between Monte-Carlo simulations 
and the semi-analytic approximation is excellent, especially throughout the 7-ray regime. 
At lower energies, our approximation remains accurate to within a factor of < 2. 



3.2. Proton Energy Losses and Equilibrium Particle Distribution 

In addition to the n^ and electron/proton-synchrotron cascades, proton synchrotron 
emission is the other ma i n contributor to the high-energy emission in hadronic models 
(JMiicke &: Protherod I2OOOI : lAharonianl I2OOOI ) . In our model, we use an asymptotic approxi- 
mation analogous to Eq. ([8]) for the proton synchrotron emissivity. 

An equilibrium proton distribution is evaluated assuming the injection of a power-law 
distribution of protons, Qpi'jp) = Qo,p7p ^'' -^(Tp! 7i,p)72,p)- As in the leptonic model de- 
scribed in ^ a self-consistent equilibrium between this injection, particle escape, and cool- 
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ing is evaluated. Unlike the case of electrons, radiative cooling time scales for protons can 
be longer than the typical dynamical time scale of the expansion of the emission region. 
We therefore account for adiabatic energy losses as well as radiative energy losses for the 
protons. 

Adiabatic losses are evaluated through 7^/7^ = —VlV^ where V is the volume of the 
emission region. Assuming a conical jet with opening angle Qj ~ l/F, the adiabatic cooling 
rate is 7p = —3 c Q^ Ip/R- The proton synchrotron cooling rate is 



Ipr 



sy 



carB^ 
Git iTipC^ 



rrir. 



Iv 



(15) 



and the photo-pion energy losses are approximated as (jAharonianll2000f ) 



Ip^py 



—cia- 



P7 



/)^ph(e*)e*7p 



(16) 



where {(Jp-yf) ~ 10^^^ cm^ is the elasticity-weighted p7 interaction cross section and e* = 
5.9 • 10^^ E^g is the energy of target photons interacting with protons of energy E at the A 
resonance, which, for pl ausible blazar param eters, is usually the most relevant contribution 
to the p7 cross section (jMiicke et al.ll2000bl ). 



The equilibrium proton distribution can then be found by integrating the analog of 
the isotropic Fokker-Planck equation (ITO|) . neglecting particle escape (which only affects the 
lowest-energy protons which do not substantially contribute to the radiative output): 



Nph 



p\ ipj 



-1 



'^Ip ^pKlp) 



(17) 



In conical jet geometries, the energy l oss rates of relativistic protons can conceivably be 
dominated by adiabatic losses (see, e.g.. lSikorall201l[ ). Because of their linear dependence on 
7p, adiabatic losses will not produce a break in the equilibrium proton spectrum. Eq. (fTTl) 
predicts a break in the proton spectrum only where radiative (proton synchroton or photo- 
pion production) losses become dominant. Figure |3] illustrates this effect. The bottom panel 
shows the proton energy loss rates due to adiabatic, synchrotron and photo-pion losses, along 
with the total cooling rate. The target photon field is dominated by the synchrotron emission 
from an electron population injected between 7e,i = 10^ and 7e,2 = 10^ with qe = 2.6, in a 
magnetic field of -B = 30 G. Protons were injected with a power-law distribution of index 
qp = 2.4, between Ep^i = 1 GeV and Ep^2 = 10^° GeV. The emission region has a radius 
of i? = 10^^ cm and expands at an angle of 6^ = l/F for F = 20. For these parameters. 
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Fig. 3. — Proton cooling rates (bottom panel) and the equilibrium proton distribution (top 
panel) evaluated according to Eq. ( IT7|) for typical synchrotron-proton-blazar model param- 
eters (see text for details). 

radiative losses (dominated by proton synchrotron) begin to dominate over adiabatic losses 
at 7p,b ~ 2 X 10^. While for all proton Lorentz factors, the cooling time scale is shorter than 
the escape time scale, a break in the equilibrium proton distribution (from index 2.4 to 3.4) 
occurs only at 7p ;,. 

The leptonic component of our hadronic jet model is being handled with the same 
procedure as described in ^ for our leptonic model. The full leptonic radiative output 
(synchrotron + SSC) is used as target photon field for photo-pion production of the hadronic 
component and for evaluating the 77 opacity of the emission region. 



3.3. Limitations 



The iKelner fc AharonianI (120081 ) templates we are using to evaluate the production 
spectra of the final decay products, neglect the synchrotron (and Comp ton) emission frorn 
intermediate decay products, i.e., muons and pions. It has been shown (JMiicke et al.ll2003l ) 
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that under certain conditions, these emissions can make a non-negligible contribution to the 
high-energy spectra in the SPB model. Neglecting these contributions is acceptable if one of 
two conditions is fulfilled: (1) The synchrotron cooling time scale of relativistic muons and 
pions are much longer than their decay time scales; or (2) proton synchrotron losses strongly 
dominate over photo-pion loses. 

For condition (1), we need to compare the decay time scales of pions (muons) with 
Lorentz factors 7,^ (7^) in the blob rest frame, 



t^a = 2.6xlO"S^s 

t', = 

(18) 



2.2 X IO^Sa^s 



to the pion (muon) synchrotron cooling time scale analogous to Eq. flT^ . This yields 
the conditions 



„ (.0 A iU *or iOI Uiont) /-,„x 

^x<! ..,,,.10^ . (19) 



7.8 X 10^^ G for pions 

5.6 X 10^° G for muons 
The Lorentz factors of secondary pions and muons are of the order of the Lorentz factors 
of the primary protons. Therefore, Eq. (fT9|) constitutes a restriction to the highest proton 
Lorentz factors allowed for a given magnetic field so that pion and muon synchrotron emission 
can be neglected. Obviously, the condition is more restrictive for muons than for pions. Our 
code automatically prints out a warning if the condition ( IT9l) is not fulfilled for either muons 
or pions. 

To check for condition (2), i.e., proton synchrotron losses being dominant over photo- 
pion losses, our code evaluates the proton synchrotron and photo-pion energy loss rates 
automatically (as plotted, for example, in Figure |3]). We only trust our results if either 
condition (1) for muons, or condition (2) is fulfilled. Specifically, this means that our semi- 
analytical hadronic model is not applicable to situations with both very high magnetic fields 
and high radiation energy densities, in which photo-pion production may dominate over 
proton synchrotron emission, and the synchrotron cooling time scales of muons (and pions) 
may be shorter than their decay time scale. 

In summary, our steady-state hadronic emission model is limited by its neglect of 
pion and muon synchrotron radiation, implying restrictions to the maximum possible field 
strength as a function of maximum proton energy injected. Protons may lose their energy 
radiatively either via the proton synchrotron channel, or photomeson production. If the 
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latter is the case, Bethe-Heitler pair production is considered subdominant, and is therefore 
neglected. The high energy photons and secondary particles from photomeson production 
may initiate cascades where the injected power is restributed by means of synchrotron radi- 
ation. We restrict ourselves to linear cascades. Inverse Compton losses are neglected, owing 
to the large magnetic-field strengths in the emission region considered here. 

External radiation fields as target field for particle and photon interactions and cascading 
are also not considered here. We point out that in the case of FSRQs, this simplification may 



not always be justified (e.g., lAtoyan fc Dermeiil2003l ). We have therefore carefully verified 



whether the external radiation fields in the 6 FSRQs modelled in this paper (see Q may 
make a substantial contribution to the target photon fields for pion production. For most of 
them, the co-moving synchrotron photon energy density in our hadronic model is at least two 
orders of magnitude larger than the co-moving-frame energy density of the external radiation 
field used for the leptonic models (which were based on the observed BLR luminosity and a 
typical size of the BLR). Direct accretion disk emission is unlikely to play a non- negligible 
role due to the unfavorable angle of incidence, which requires very high proton Lorentz 
factors for photo-pion production, as would the interaction with an infrared radiation field 
from a dust torus. However, for 3C454.3 and PKS 0528-1-134, we find that the BLR radiation 
field (Doppler-boosted into the co-moving frame) may be of the order of or even larger than 
the synchrotron photon energy density, depending on the (unknown) radius of the BLR and 
the location of the 7-ray emitting region with respect to the BLR. We therefore caution that 
our neglect of external radiation fields may not be valid in those two cases. 



4. Comparative Modeling of Fermi-LAT detected blazars 

In this section, we will apply our leptonic and hadronic models descr i bed in the previous 



sections, to a set of 7-ray blazars detected by Fermi-LAT . lAbdo et al.l ((20101) presented a 
comprehensive compilation of simultaneous or contemporaneous broadband SEDs of a large 
sample of Fermi-LAT detected blazars. Out of their list, we selected a subset of blazars which 
fulfilled the following criteria: (1) good simultaneous broadband SED coverage, including at 
least optical/IR, X-ray, and 7-ray data; (2) known redshift; (3) information about short-term 
variability; (4) apparent superluminal motion; (5) in the case of FSRQs: information about 
the accretion-disk and BLR luminosity. Condition (1) is necessary to allow for meaningful 
SED fitting. Condition (2) is required to determine the overall luminosity requirements of 
models. Condition (3) allows us to derive constraints on the size of the emission region 
from the observed variability time scale tvar via R < Dct^a_r/{^ + z); (4) superluminal 
motion with a speed /3±,app allows us to place a lower limit on the bulk Lorentz factor of 
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r > \//3iapp + l5 ^^'^ (^) knowledge of the accretion- disk and BLR luminosities further 
constrains the modeling of FSRQs in which these sources of external radiation fields are 
generally believed to be important. 

We found that the following objects fulfilled our criteria: The FSRQs 3C 273, 3C 279, 
3C 454.3, PKS 1510-089, PKS 0420-01, and PKS 0528+134; the LBLs BL Lacertae, AO 0235+164, 
S5 0716+714, and OJ 287; and the IB Ls W Comae and 3C 66A. SED data for these objects 



are obtained from lAbdo et al.l (120101 ). and Table [T] lists observational data that are used to 



further constrain our models. 

Our fitting procedure is a "fit-by-eye" method, starting with plausible values for the pa- 
rameters that are not directly constrained by observations. The unconstrained parameters 
are then adjusted to obtain an acceptable representation of the SED. Simpler models with 
fewer parameters, such a s a single- zone, static SSC model, allow for a detailed x^ minimiza- 



tion technique (see, e.g.. iMankuzhiyil et al.l 120111 ). However, given the substantial number 



of adjustable parameters, such a fitting method is infeasible for the models considered here. 
While we are confident that our best-fit parameters are in the correct ball-park range for 
these objects, the lack of a rigorous x^ minimization strategy prevents us from determining 
errors on the fit parameters. Therefore, while we will proceed to a global assessment of 
parameter values among different classes of blazars, we will refrain from any quantitative 
statements in this respect. 

In the course of the fitting, we strive to achieve acceptable fits with parameters close to 
equipartition between the dominant particle species (electrons in the leptonic model, protons 
in the hadronic model) and the magnetic field. This is motivated by two arguments: (1) 
If the relativistic jets ofAGN are powered by rotational energy from the central black hole 



( iBlandford fc Znajekl 119771 ) . the jets are expected to be initially Poynting-fiux dominated, 
and the energy carried in electromagnetic fields needs to be transferred to relativistic particles 
in order to produce the observed high-energy emission. This energy conversion is expected 
to cease as approximate equipartition between the magnetic field and relativistic particles 
is reached, so that the jets are not expected to become matter dominated in the central few 
parsecs of the AG N, where the high-energy emission in blazars is believed to be produced 



( jLyubarskyl |2010| ) . (2) If magnetic pressure plays an essential role in collimating AGN jets 
out to kpc scales, the particle pressure can not largely dominate over the magnetic pressure 
in the inner few parsecs of the AGN. For these reasons, we prefer model parameters with 
magnetic fields dominating the pressure and energy density in the emission region over 
particle-dominated scenarios. 

For all objects, we produced one leptonic and one hadronic model fit to the contem- 
poraneous SED. In the case of 3C 66A, the catalog value of the redshift oi z = 0.444 is 
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Object Name Type 



/5j 



.,app 



f-var 



^disk [erg s ^] Lblr [erg s ^] 



3C 273 


FSRQ 


0.158 


13(1) 


Id (2) 


1.3 X 10^7 (3) 


9.1 X 10^5 (4) 


3C 279 


FSRQ 


0.536 


20.6 (5) 


2d (6) 


2.0 X 10^5 (7) 


2.0 X 10^4 (8) 


3C 454.3 


FSRQ 


0.859 


14.9 (^) 


1 hr (9) 


1.7 X 10^7 (10) 


2.5 X 10^5 (11) 


PKS 1510-089 


FSRQ 


0.36 


20.3 (5) 


1 d n 


1.0 X 10^6 (13) 


— 


PKS 0420-01 


FSRQ 


0.916 


7.5 (5) 


4 d n 


1.5 X 10^6 (11) 


4.3 X 10^^ (11) 


PKS 0528+134 


FSRQ 


2.06 


30 (15) 


1 d n 


1.7 X 10^7 (1*^) 


1.8 X 10^6 (16) 


BL Lacertae 


LBL 


0.069 


10.7 (5) 


1.5 hr (17) 


6.0 X 10^4 (18) 


6.3 X 10^2 (11) 


AO 0235+164 


LBL 


0.94 


2(5) 


1.5 d (19) 


3.4 X 10^4 (11) 


3.7 X 10^3 (11) 


S5 0716+714 


LBL 


0.31 


10.2 (5) 


15 mill (20) 


— 


— 


OJ287 


LBL 


0.306 


15.3 (5) 


2 hr n 


1.1 X 10^5 (11) 


2.7 X 10^3 (11) 



5 hr (23) 

6 hr (25) 



W Comae 
3C 66A 



IBL 
IBL 



0.102 
0.44 ? 



2(22) 

27 (24) 



7.2 X 10^4 (11) 5.0 X 10^1 (11) 



Table 1: Observational data used to constrain ou r mod els. S uperscripts in paranethe ses 
re fer to the following; referenc es: 1 



Bottcher et al. 



J2008bl): 10 



feoOTl): 



Raiteri et al. 



Pucella et al. 



Palma et al. 



Lister et al.l ( 120091): 2 = ICourvoisier et a. 



Vasudevan fc FabianI (120091): 4 = IPeterson et al.l (120041: 5 = iHovatta et al. (120091): 6 



Hartman et al. 



J2008ch: 11 



(120011): 8 



Xie et a 



teOOSi); 14 = iD'Arcaneelo et al.l (120071): 15 



(2011): 17=|Vi 



(l2008al): 20 =ISasada et al.l (120081): 21 =lFan et al.l (120091): 22 



lataet al.l (2002): 18 



Plan et al 



JiopsI); 



Tagliaferri et al.l tOOOi ): 24 = I.Torstad et al.l ( 120011 ): 25 = iTakalo et al. 



12 



1119881); 3 



(120051): 9 =lRaiteriet al.l 



Marscher et al. 



Raiteri et a. 



J2OI0I); 



orstad et al.l (120051): 16 



13 



(120091): 19 = Raiteri et al. 



Massaroet al.l (I2OOII): 23 



(Il99d ) 
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Fig. 4. — Leptonic model fits to the 6 FSRQs in our sample. See table |2] for parameters. 
Dotted = synchrotron; dashed = accretion disk; dot-dashed = SSC; dot-dash-dashed = EC 
(disk); dot-dot-dashed = EC (BLR). 
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highly questionable. The recent analyses of lPrandini et al.l (120101) and I Ab do et al.l (l201l[ ) 
place the object at a likely redshift of z ~ 0.2 - 0.3, while iFurniss et al.l ( 120131 ) provide a 
lower limit on the redshift of 2; > 0.3347, based on UV absorption features by intervening 
intergalactic medium. In our modeling, we assume a fiducial redshift oi z = 0.3 for 3C 66A. 
While, of course, the parameter val ues slightly change w hen using a slightly higher redshift 
(say, z = 0.35 to be consistent with iFurniss et al.l ( l2013l )). the overall conclusions from the 
modeling remain unaffected. All our model SED s are correc t ed for 77 absorption by the 
extragalactic background light using the model of iFinke et al.l (I2OIOI ). 



4.1. Leptonic Model Fits 

Figures H] - [6] show the SEDs and our leptonic model fits for the 12 blazars we selected 
for this study. The fit parameters used and a few quantities derived from those parameters, 
are listed in Table HI In general, satisfactory fits to most blazar SEDs with parameters close 
to (within an order of magnitude of) equipartition can be achieved with the leptonic model 
described above. 

In agreement with many previous studies, we find that FSRQs require a dominant 
contribution from external-Compton (on the direct accretion disk and the BLR emission) 



in order to provide acceptable SED fits (e. g., ISambruna et al.l Il997l : iGhisellini et al.l Il998 



Mukherjee et al.lll999l : iHartman et al.ll200ll ). In our model fits, we are able to produce the 
spectral breaks observed in the Fermi-LAT spectra of many blazars, with a superposition 
of two different 7-ray emis sion components, as previously proposed for the case of 3C454.3 
by iFinke fc Dermerl (120101 ). Typical minimum variability time scales from light-travel time 
arguments are of the order of 1 - 2 days, consistent with the day-scale variability seen in 
most Fermi-detected FSRQs. 



As argued previously, e.g. by iMadejski et al.l (Il999[ ) and iBottcher &: Bloom! (I2OOOI ) for 
the case of BL Lacertae, also LBLs are represented with more plausible parameters when 
including an EC component using a low-luminosity accretion disk and BLR, compared to a 
pure SSC model. In the case of S5 0716-1-714, where no observationally motivated estimates 
of the luminosity of the accretion disk or the BLR could be found, we have assumed the 
existence of a dust torus with a very low infrared luminosity, consistent with the absence 
of any direct observational evidence for it. Our models allow for variability on time scales 
of a few hours for the LBLs modelled here. While the SEDs of BL Lacertae, S5 0716+714 
and OJ 287 are well represented by a leptonic SSC+EC model, our model has problems 
representing the very hard 7-ray spectrum of AO 0235+164 above a few GeV as our model 
7-ray spectra are truncated by Klein-Nishina effects. 



Object 


7e,l 


7e,2 


Qe 


B[G] 


D 


/blr 


Li 




r^ 

^B 


ese 


var,mm 


3C273 


1 X 10^ 


5 X 10^ 


3.5 


2.0 


13 


5.4 X 10-^ 


6.0 


130 


0.63 


0.11 


4.1 


3C279 


1 X 10^ 


1 X 10^ 


3.0 


0.7 


17 


1.7 X 10-2 


5.8 


29 


5.4 


0.94 


29 


3C 454.3 


800 


5 X 10^ 


3.0 


2.1 


15 


8.0 X 10-2 


22 


870 


75 


3.5 


52 


PKS 1510-089 


1 X 10^ 


2 X 10^ 


3.1 


0.8 


20 


4.1 X 10-3 


2.1 


41 


2.1 


1.0 


9.4 


PKS 0420-01 


1.4 X 10^ 


5 X 10^ 


3.4 


2.5 


8.0 


1.0 X 10-^ 


6.9 


670 


38 


5.4 


110 


PKS 0528-M34 


700 


1 X 10^ 


3.0 


3.0 


19 


9.1 X 10-3 


9.4 


600 


110 


11 


45 


BL Lacertae 


1.1 X 10^ 


1 X 10^ 


3.2 


2.5 


15 


1.3 X 10-1 


0.44 


4.6 


0.41 


0.94 


1.8 


AO 0235+164 


700 


5 X 10^ 


3.7 


1.3 


25 


1.1 


12 


87 


33 


2.8 


21 


S5 0716+714 


1.9 X 10^ 


5.5 X lO'' 


3.8 


3.8 


20 


d 


1.3 


6.7 


14 


10 


4.9 


OJ287 


800 


5 X 10^ 


3.8 


3.5 


15 


1.5 X 10-1 


1.8 


16 


15 


8.2 


9.7 


W Comae 


1 X 10^ 


8 X 10^ 


2.4 


1.5 


30 


2.8 X 10-2 


0.30 


1.52 


0.30 


1.0 


0.68 


3C 66A 


9 X 10^ 


3 X 10^ 


2.8 


0.065 


40 


d 


13 


8.2 


0.45 


0.034 


13 



Table 2: Parameters for the leptonic SED fits shown in Figures H] - [6l " The factor /blr, 
determining the energy density of the re-processed disk radiation field, is defined as /blr = 
'7blr/-Rblrj values are in pc"^. ^ Powers are in units of lO^"^ erg s~^; '^ minimum allowed 
variability time scale in hours. '^ For S5 0716+714 and 3C 66A, no accretion disk luminosity 
has been measured/constrained; For S5 0716+714, an isotropic external radiation field with 
Mext = 5 X 10"^ erg cm"^ and Tbb = 2 x 10'' K has been used for the SED fit; for 3C 66A: 
Mext = 1-3 X 10-8 erg cm'^ and Tbb = 10^ K. 
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Fig. 5. — Leptonic model fits to the 4 LBLs in our sample. See table |2] for parameters. 
Dotted = synchrotron; dashed = accretion disk; dot-dashed = SSC; dot-dash-dashed = EC 
(disk); dot-dot-dashed = EC (BLR). 
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In blazars in which we find that a substantial contribution from EC on direct accretion- 
disk emission is relevant, we find a best-fit distance of the 7-ray emission region from the 
central black hole in the range zq < 0.1 pc. An important contribution from EC on BLR 
emission would imply a characteristic distance of z < 1 pc, while a much larger distance 
would be consistent with situations in which EC on infrared emission (from warm dust) 
dominates the 7-ray production. As described in ^ the external radiation fields from the 
BLR and/or warm dust are modelled as isotropic in the AGN rest frame so that for such 
fits, the distance to the black hole is not necessary input parameter for our model. 

While BL Lac objects detected at > 100 GeV 7-rays by ground-based Atmospheric 
Cherenkov Telescopes have often been found to be well represented by pure SSC models, 
the careful ai ialysis of the simult a neous SEDs of th e VERITAS-detected IBLs W Comae 



and 3C 66A (jAcciari et al.ll2009bl : lAbdo et al.ll201ll ) revealed that also in these two cases 



a pure SSC fit would require rather extreme parameters with magnetic fields far below 
equipartition. This situation could be remedied allowing for a contribution from external 
Compton to the 7-ray emission. These findings are confirmed in this study. However, the 
unusual, apparent upward curvature of the Fermi-LAT spectrum of 3C 66A and W Comae 
can not be satisfactorily represented with our model, irrespective of the dominant target 
photon field for external Compton scattering. 

We find that the two IBLs in our sample require systematically higher Doppler factors 
than the LBLs and FSRQs. Also, as expected, BL Lac objects are characterized by less 
powerful jets (i.e., lower Lg) than FSRQs, and the two IBLs require harder electron spectra 
than the LBLs and FSRQs. We do not find a systematic difference in magnetic-field values 
between BL Lacs and FSRQs, neither are the characteristic electron energies (represented 
by 7e,i) systematically different between the two classes of objects. If charge neutrality in 
blazar jets is provided by cold protons (rather than positrons), our fits indicate that the 
kinetic energy carried by the jets should be dominated by protons. 



4.2. Hadronic Model Fits 

Overall, we find that the SEDs of BL Lac objects, both IBLs (Fig. E) and LBLs (Fig. [8]), 
can be well represented with our model. The typically low ratio of 7-ray to X-ray fiux and 
the hard 7-ray spectra observed in BL Lac objects (as compared to FSRQs) is naturally 
obtained by photo-pion induced cascade emission with substantial contributions from proton- 
synchrotron radiation. In some cases, an additional contribution from the leptonic SSC 
emission aids in producing the observed X-ray flux. However, even though the addition 
of photopion induced cascade emission to a proton synchrotron component is, in principle. 
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0.27 




1.2 


X 10-* 


9.7 


PKS 0528-f 134 


150 


1.0 X lO" 


3.8 


30 


20 


1.1 


X 


10» 


2.0 


1.9 


44 


1020 


4.4 X 10" 


-:i 


4.3 


X 10-6 


17 


BL Laccrtac 


700 


1.5 X 10* 


3.5 


10 


15 


1.9 


X 


lo" 


2.4 


0.087 


9.8 


39 


3.4 X 10" 




8.9 


X 10-' 


1.3 


0235+164 


200 


750 


3.0 


15 


25 


4.3 


X 


10" 


1.9 


1.5 


10 


130 


1.9 X 10" 


-3 


1.5 


X 10-^ 


4.3 


S5 0716-1-714 


900 


3.0 X 10* 


2.9 


20 


15 


2.7 


X 


10" 


2.0 


0.089 


0.14 


1.4 X 10* 


0.85 




6.2 


X 10-^ 


15 


O.T 287 


350 


4.0 X 10* 


4.1 


20 


15 


1.0 


X 


10" 


1.6 


0.53 


8.3 


66 


0.042 




6.3 


X 10-* 


2.6 


W Comae 


800 


2.1 X 10* 


2.6 


30 


15 


1.9 


X 


lo" 


2.0 


0.014 


0.021 


560 


0.037 




6.6 


X 10-=' 


0.68 


3C 66A 


750 


1.3 X 10* 


2.8 


10 


30 


1.2 


X 


10" 


2.0 


0.32 


1.2 


24 


6.5 X 10" 


-4 


2.7 


X 10-^ 


0.57 



Table 3: Parameters for the hadronic SED fits shown in Figures [7] - [91 " Magnetic field in 
units of Gauss. '' Kinetic luminosity in relativistic electrons in units of 10''^ erg s~^. '^ Kinetic 
luminosity in relativistic protons in units of 10^* erg s~^. '^ Partition fractions defined as 
eij = Li/ Lj. ^ Minimum allowed variability time scale in hours. 
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able to reproduce concave 7-ray spectra, the steep upturn of the Fermi-LAT spectrum of 
3C 66A and W Comae at GeV energies is still not well represented also by this model, 
and may indicate the need for an additional emission component (possibly from muon/pion 
synchrotron radiation, with consequences for the values of the fit parameters). We note that 
the required power in relativistic protons, Lp, is very large, in the range ~ lO'^^ - lO^^erg 
s~^, which is significantly higher than the observed radiative luminosities of these objects. 

Fits within the framework of the hadronic model to the SEDs of the FSRQs in our 
sample are shown in Fig. [71 Although the hadronic model appears to satisfactorally fit 
most of the snap-shot SEDs, the typically observed spectral break at GeV energies seems 
problematic in the cases of 3C 273 and 3C 279. 

Also here, the 7-ray component is dominated by proton synchrotron radiation, with 
some contribution from proton initiated pair cascade emission above ~ 10 GeV. The too 
steep decline of the proton synchrotron emission at a few GeV may indicate too steep a 
decline of the radiating proton distribution in the cutoff region (here considered a hard- 
cutoff), and hence possibly an im pact of the acceleration mechanism in this region (e.g.. 



Protherod I2004J : lAharonianl I2OOOI ) . The observed hard X-ray to soft 7-ray spectral shapes 
can be fitted if proton synchrotron emission is the dominant radiation channel (in order not 
to over-predict the X-ray flux due to reprocessing of UHE 7-ray emission in cascades). This 
indicates the presence of high magnetic field strengths in the emission region. 

The hadronic FSRQ fits also require extreme jet powers, in some cases exceeding lO^^erg 
s~^ in relativistic protons. 

In summary, we find that the hadronic model presented here provides satisfactory fits to 
most of the bright blazar SEDs of our sample. However, the declining (in uF^) broken power- 
law shape observed in many Fermi-LAT spectra of FSRQs causes problems for adequate 
fits in two cases within the limits of this model. Fits to most objects have been achieved 
with magnetic fields in the range i? ~ 10 - 30 G. All fits required proton acceleration to 
energies of Ep > 10^^~^^ eV. Nearly all model fits used a strong cascade component to aid 
modeling the high energy GeV data. This leads to the requirement of a disproportionately 
large proton luminosity as compared to the magnetic field luminosity. As a consequence in 
all hadronic model fit parameter sets presented here the total jet power is dominated by 
protons. Alte rnatively, hadronic models that take charged 7r//i synchrotron radiation into 



account (e.g., iMiicke et al.l 120031 ) could lower the strength of the pair casade component, 
thereby reducing the proton contribution to the jet power at the expense of magnetic field 
power. In order to produce the observed (electron-synchrotron) IR - UV emission, the model 
requires lower characteristic electron energies for FSRQs (7e,i ~ 100) than for BL Lac objects 
(7e,i ^ 10^), along with harder electron spectra in BL Lac objects. All our fits were achieved 
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with characteristic Doppler factors of D ~ 10 - 30. 



5. Summary and Conclusions 

In this paper, we have described the development of new implementations of stationary, 
single-zone leptonic and hadronic models. Our leptonic model allows for arbitrary external 
photon sources, and solves self-consistently for an equilibrium between relativistic particle 
acceleration (iniection), rad iative cooling, and escape. Our hadronic model is based on the 



Kelner fc Aharoniad (120081 ) templates for the final products of photo-pion production, and 
uses a new, semi-analytical method for calculating the output from ultra-high-energy induced 
pair cascades. 

We have used both leptonic and hadronic models to fit the contemporaneous SEDs of 12 
Fermi-LAT deteced blazars with good multiwavelength coverage and additional observational 
constraints on model parameters. We find that the SEDs of all types of blazars can be 
well represented with leptonic models with parameters close to equipartition between the 
magnetic field and relativistic electrons in the emission region. However, our leptonic model is 
unable to provide a good fit to the hard Fermi-LAT spectrum of AO 0135+164. The problem 
lies in the mismatch between the very steep synchrotron (IR - optical - UV) continuum, 
as opposed to the very hard 7-ray spectrum, and Klein-Nishina effects at the highest 7- 
ray energies. We confirm that even intermediate BL Lac objects are more appropriately 
fit including an external radiation field as source for Compton scattering to produce the 
observed 7-ray emission. FSRQs are characterized by systematically more powerful jets, 
but lower bulk Lorentz factor and softer electron spectra than BL Lac objects. If charge 
neutrality in blazar jets is provided by cold protons (rather then relativistic positrons), our 
fits indicate that those protons are dominating the kinetic power carried by the jets by about 
an order of magnitude. 

The hadronic model presented here has difficulty describing the GeV-break in the SEDs 
of two FSRQs, but provides appropriate fits for all other blazars in our sample. However, 
the fits require very large powers in relativistic protons, of Lp ~ 10^^ - lO'^^ erg s~^, in most 
cases dominating the total power in the jet. 

We thank Paolo Giommi for sending us the SED data for our modeling, and the anony- 
mous referee for a helpful and constructive report. MB acknowledges support by the South 
African Department of Science and Technology through the National Research Foundation 
under NRF SARChI Chair grant no. 64789. This work was supported by NASA through 
Astrophysics Theory Program Award NNX10AC79G and Fermi Guest Investigator Grant 
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W Comae 



3C66A 





Fig. 6. — Leptonic model fits to the 2 IBLs in our sample. See table |2] for parameters. 
Dotted = synchrotron; dashed = accretion disk; dot-dashed = SSC; dot-dash-dashed = EC 
(disk); dot-dot-dashed = EC (BLR). 
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Fig. 7. — Hadronic model fits to the 6 FSRQs in our sample. See table [3] for parameters. 
Dotted = electron-synchrotron; dashed = accretion disk; dot-dashed = SSC; dot-dot-dashed 
= proton-synchrotron. 
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Fig. 8. — Hadronic model fits to the 4 LBLs in our sample. See table |3] for parameters. 
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Fig. 9. — Hadronic model fits to the 2 IBLs in our sample. See table |3] for parameters. 



